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Abstract. A detection of the stacked integrated Sachs- Wolfe (ISW) signal in the CMB of rare 
superstructures identified in the SDSS Luminous Red Galaxy catalogue has been reported 
at very high statistical significance. The magnitude of the observed signal has previously 
been argued to be more than 3a larger than the theoretical ACDM expectation. However, 
this calculation was made in the linear approximation, and relied on assumptions that may 
potentially have caused the ACDM expectation to be underestimated. Here we update the 
theoretical model calculation and compare it with an analysis of ISW maps obtained from 
A-body simulations of a ACDM universe. The differences between model predictions and 
the map analyses are found to be small and cannot explain the discrepancy with observation, 
which remains at > 3er significance. We discuss the cosmological significance of this anomaly 
and speculate on the potential of alternative models to explain it. 
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1 Introduction 

The late-time integrated Sachs- Wolfe (ISW) effect can be used as a cosmological probe that is 
sensitive to the dynamical effects of dark energy [1]. It is manifested as secondary anisotropies 
in the cosmic microwave background (CMB) radiation which are introduced as photons from 
the last scattering surface travel through time-evolving gravitational potentials [2]. Under 
the assumption of spatial flatness and at linear order, the late-time evolution of potentials 
only occurs in an accelerating Universe. The detection of the ISW effect at linear order is 
therefore an important independent confirmation of the effects of dark energy, and can be 
used to test the standard A Cold Dark Matter (ACDM) cosmological model. 

Unfortunately, the amplitude of the ISW effect is an order of magnitude smaller than 
the intrinsic fluctuations in the CMB. Therefore detection of the signal requires careful cross- 
correlation of the CMB with data from large-scale structure catalogues that trace the matter 
density distribution [3, 4]. Several such studies have been performed using different tracers, 
e.g. [5-17]. A range of different results is obtained; for a summary, see [18]. The use of a com- 
bination of several large-scale structure data sets has been found to give overall significance 
of ~ 4cr for the observation of a cross-correlation [19-21]. Greater sensitivity is required to 
test whether the form of the cross power spectrum is as expected in ACDM, though there 
are intriguing hints that the amplitude of cross-correlation with the SDSS LRG catalogue 
may be larger than expected [19-22], as well from the WISE catalogue [23]. 

An alternative approach to the detection of the ISW signal has been to focus on the 
contributions from individual superstructures alone [24, 25]. In this approach, the authors 
identified the most extreme over- and under-dense superstructures in distribution of LRGs 
in the SDSS DR6 MegaZ photometric catalogue, and then measured the CMB temperature 
along the lines of sight corresponding to the directions of the identified structures. By stacking 
the CMB images so obtained and filtering with a compensated top-hat filter, they reported 
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detections of hot spots correlated with overdensities and cold spots with underdensities, 
which may be attributed to the ISW effect. The amplitude of the correlation was found to be 
significant at > 3a each for hot and cold spots individually, and at ~ 4.4<r for both combined. 
This is the highest significance detection yet obtained using a single tracer population. 

However, the size of this detected signal is difficult to understand within the ACDM 
model [26, 27]. A recent analysis in ref. [28] showed that the reported stacked signal is 
more than 3 standard deviations larger than the maximum expectation in the ACDM model 
if the primordial density perturbations were Gaussian, even with optimistic assumptions. 
The analytic calculation in [28] was however fully linear 1 and based on certain theoretical 
assumptions including the spherical nature of the average density fluctuation, and therefore 
did not take into account a potential biasing effect arising from smearing of the LRG sample in 
the radial direction due to photometric redshift errors. The effects of these simplifications are 
described in greater detail in section 3.2. These simplifications might in principle have affected 
the predicted maximum value of the stacked signal. This left open the small possibility 
that the most optimistic prediction from the ACDM model might still be compatible with 
observation. 

To test the quantitative effect of these simplifications on the calculation in ref. [28] 
the predictions of the analytic model should be compared with results obtained from large 
A-body simulations. To our knowledge, there does not exist an A^-body simulation of the 
ISW effect in ACDM which covers the same redshift range and angular footprint as the 
SDSS photometric survey. However, ref. [29] provides two realizations of sky maps of the 
ISW temperature anisotropics, including non- linear contributions, calculated for A^-body 
simulations based on the ACDM model of the matter distribution in two different redshift 
intervals, corresponding to two slabs of depth 1 h~ 1 Gpc. By filtering and stacking images 
from these maps in the same manner as done originally in [24, 25] we estimate the size 
of the expected maximum possible signal in ACDM without recourse to the assumptions 
made in the previous analytic treatment. Such a treatment gives an upper bound on the 
maximum possible ISW signal from superstructures and is unlikely to be achieved in any 
realistic observation. 

Comparison with this upper bound obtained from the simulated maps shows that the 
spherical model slightly overestimates the maximum possible signal at large redshifts and 
underestimates it at small redshifts. The numerical difference between the spherical model 
prediction and the theoretical maximum is however always small. In particular, it is not 
sufficient to explain the discrepancy between the observed stacked temperature signal and 
the ACDM expectation obtained from the spherical model. 

Using the simulated ISW maps, we are also able to investigate the variation of the 
maximum possible signal-to-noise ratio with the number A" of stacked structures used and 
the radius Of of the applied compensated top-hat filter. We find that the stacked ISW 
signal of superstructures in a ACDM universe with Gaussian primordial perturbations is not 
observable for any choice of A" and Of, unless the survey volume is much larger than the 
SDSS survey. We conclude that the high significance detection reported in ref. [24] is not 
compatible with the ISW effect in such a cosmology, and the conclusions of ref. [28] are 
confirmed. 

The structure of the paper is as follows. In section 2 we briefly introduce the necessary 
background theory and describe the methodology used in the observation of the stacked 

1 Ref. [27] accounted for some non-linear effects but used a less conservative method for calculating the 
expectation value of the signal. The qualitative conclusions were the same. 
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signal. We then turn to the analytic spherical model. The main features of the model 
and the improvements over ref. [28] are described in section 3.1, with several mathematical 
details relegated to Appendix A. We explain the assumptions made and potential limitations 
of the model in section 3.2. Then in section 4 we introduce the ISW maps produced from 
iV-body simulations and describe our method of analysis. The results of this analysis and 
a comparison with the spherical model are presented in section 5. Finally we conclude the 
paper with a discussion of the implications of the result in section 6. 

2 The ISW signal of superstructures 
2.1 Theory 

The ISW effect is a secondary temperature anisotropy that arises when photons from the 
CMB travel through time-evolving gravitational potentials. At late times in the standard 
flat cosmology, this time evolution arises due to the effect of the non-zero dark energy density 
The temperature anisotropy AT(n) thus induced is [2] 
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AT(n) = -^T I <5>{r,z,h)adr, (2.1) 
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where To is the mean CMB temperature, tl is the radial comoving distance to the last 
scattering surface, a is the scale factor, <1> is the time derivative of the gravitational potential 
and c is the speed of light. 

The potential # is related to the density fluctuation in the conformal Newtonian gauge 
via the Poisson equation, 

V 2 $(x, t) = 4vrGp(t)a 2 (5(x, t), (2.2) 

where p(t) is the mean matter density of the universe and 5 = (p — p)/p is the density 
contrast. Changing to Fourier space and taking the time derivative gives 
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In ref. [29], the continuity equation 5(k,t)+ik-p(k,t) = is used in eq. (2.3) to produce 
the simulated ISW maps, where p(k,t) is the Fourier transform of the momentum density 
divided by the mean mass density, and can be obtained directly from the simulation. The 
analytic model of ref. [28] and its slightly modified version introduced in section 3 both use 
the linear growth approximation 5(k, t) = D(t)5(\i,z = 0), where D(t) is the linear growth 
factor. Substituting this into eq. (2.3) yields: 

*(k,z) = \ (^j 2 n m H ^[l-(3(z)]5(k,z), (2.4) 

where (3(t) = din D/ din a denotes the linear growth rate. As this model also assumes spher- 
ical symmetry of the structures, eq. (2.4) may be rewritten in real space as [28] 

<f>(r,z) = ~n m H$G(z)F(r), (2.5) 
where G(z) = H(z)(l — f3{z))D{z) / a is the ISW linear growth factor and 

rr 12 i-oo 

F{r) = / — S(r')dr' + / r'6(r')dr' , (2.6) 

JO r Jr 

with S(r') evaluated at z = 0. 
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2.2 Observation of the stacked signal 

The signal reported in ref. [24] was obtained in the following way. The authors used a 
sample of 1.1 million luminous red galaxies (LRGs) with photometric redshifts in the range 
0.4 < z < 0.75, and a median redshift of z ~ 0.52, from the SDSS release DR6 [30], covering 
a footprint of 7500 square degrees about the North Galactic pole. Within this sample, they 
identified over- and under-dense fluctuations in the galactic density field using the structure- 
finding algorithms VOBOZ [31] and ZOBOV [32] respectively. The structures thus found were 
ranked according to their density contrast; this is roughly correlated with both the likelihood 
of such a structure existing and the relative ISW signal produced by it. 

Given such a ranking, the N top-ranked "superclusters" and "supervoids" in the sample 
were selected, and the CMB images in the directions along the lines of sight of these struc- 
tures were stacked and averaged. The CMB dataset used was an inverse-variance weighted 
combination of the WMAP 5-year Q, V and W maps [33] with Galactic foreground template 
maps subtracted and the KQ75 mask applied. The temperature in the stacked images was 
then averaged with a circular compensated top-hat filter of the form 



where ou t = V20f- The effect of the filter is thus to take the average temperature within 
an inner circle of radius Of and subtract from it the average temperature within an outer 
circular ring of equal area, in order to remove the effect of primary CMB fluctuations on 
scales larger than the filter width. The noise due to the primary CMB fluctuations was 
estimated by randomising the positions of the stacked directions on the actual CMB map, 
and by using model CMB maps smoothed to the WMAP resolution. 

For the choice N = 50 and 9f = 4°, the result obtained was AT C = 7.9 db 3.1 /xK for the 
hot spots associated with superclusters, AT V = — 11.3 ± 3.1 [aK for the cold spots associated 
with supervoids, and AT = 9.6 ± 2.2 fiK for the two combined (clusters minus voids), i.e., a 
4.4<7 significant detection [24]. The variation in signal strength with change in N and Of was 
also investigated and found to be maximum at N = 50 and Of = 4° . It is noteworthy that 
the signal obtained is not frequency dependent [24] and is therefore unlikely to be caused by 
contamination by foregrounds or by extragalactic radio sources. 

A follow-up investigation by the same authors [22] using a template-fit analysis with 
an ISW template reconstructed from the LRG density field confirmed the existence of the 
high-significance correlation between the observed supercluster and supervoid locations and 
the hot and cold spots on the CMB. This study also did not find any evidence of foreground 
contamination. 

3 The analytic spherical model 

3.1 Calculation of the expected signal 

The analytic spherical model for the stacked ISW signal of superstructures introduced in 
ref. [28] is based on the statistics of peaks of a Gaussian-distributed random field. We 
assume that the matter density contrast 6, smoothed on a length-scale Rf large enough that 
all fluctuations are within the linear regime, is such a random field. Then by using the 
corresponding matter power spectrum for the ACDM cosmological model, one can calculate 




(2.7) 
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the predicted (comoving) number densities of different populations of maxima (and minima) 
in S. The differential number density of maxima is [34] 

r , v e^ 2 /2 exp [_( x _ 7I/ )2/2(l -7 2 )1 . 
A/" max (^)dz/dx = [27r(1 _ 72)] i/2 dl/dx ' (3 - 1} 

where f = 5/<7o> x = —V 2 S/a2, Co and 02 are members of a set of spectral moments of the 
power spectrum of the density field, and 7 and R* are spectral parameters related to these 
moments. Definitions for these quantities as well as the full closed-form expression for f(x) 
are provided in Appendix A. 

We then identify the points of maxima in S with the overdense structures in the galaxy 
density field found in the SDSS LRG catalogue and points of minima with voids. The galaxy 
density field S g provides a biased tracer of 5. This could be modelled using a bias relationship 
6g = bS, but in order to obtain the most optimistic estimate of the maximum achievable ISW 
signal it is sufficient to assume simply that the largest fluctuations in S g occur at the same 
locations as the largest fluctuations in S. 

If a maximum (minimum) of given height (depth) So = voo and with a given central 
x exists at a point r = 0, the mean profile of the density contrast S(r; v, x) about that 
point can be calculated [34]. In general the structures about a point of extremum have a 
triaxial ellipsoid shape, but in the absence of any preferred direction one can average over all 
orientations of the axes and obtain a mean spherical profile about a peak [34, 35]: 



1 f°° k 2 sm(kr 



X7 (x — ^v)R 2 k 2 



+ 



dk, (3.2) 



(1 _ 7 2) 3 7 (1 -7 2 ) 

where P{(k) denotes the matter power spectrum filtered on the scale Rf using a spherical 
Gaussian filter. As long as fluctuations remain within the linear regime, the symmetry 
between peaks and voids is maintained. This can be seen from the fact that eqs. (3.1) and 
(3.2) are unchanged under a change of sign of v and x for voids. 

We model the structures identified in the SDSS LRG catalogue as being due to spheri- 
cally symmetric peaks or voids in the matter density field, with profiles of the form eq. (3.2). 
For a given choice of £0 and x we calculate the ISW temperature shift AT(6;So,x) caused 
by a structure with such a profile by placing the centre of the sphere at a redshift z = 0.52 
(corresponding to the median redshift of the LRG catalogue) and using eqs. (2.5) and (2.6) to 
perform the ISW line-of-sight integral of eq. (2.1) for any photon incidence angle through 
the structure. To this temperature profile we then apply the same compensated top-hat filter 
of width Op (eq. (2.7)) used in [24]. The resulting quantity AT (So, x; Of) is the observationally 
relevant average temperature signal due to this structure. 

The evaluation of J\f max (v, %) and AT (So, x; Of) for the ACDM model requires the spec- 
ification of the cosmological parameters. We use the same parameter values as in ref. [28], 
which are the mean values obtained from a fit to WMAP 7-year [36] and SDSS DR7 data [37]. 
However, the precise choice is unimportant, since variation of parameter values within the 
ranges allowed by the concordance cosmology has a negligible effect on the spherical model 
prediction [28]. To calculate the matter power spectrum P(k) at the requisite redshift we 
use CAMB [38]. 

To obtain an expectation value of the average signal caused by a sample of structures of 
different central density contrasts and radii 2 we now take a weighted average of AT (So, x; 9p) 



2 Note that the radius of a structure with given So is determined by the normalized second derivative of 
the density field at the centre, x = — V 2 5/o"2- Larger values of x correspond to narrower profiles. 
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over the range of values of Sq and x, with the weighting function chosen according to the 
relative number of such structures present in the sample of structures. Thus for a sample of 
the N largest and most overdense structures contained within some comoving volume V, we 
obtain the expectation value 

v r 1 r Xcut 

(AT pk (N,e F )} = - / ATp k (5 ,x;^)Mnaxao d<J odx, (3.3) 

iv y<5§ Jo 

where 5q and x cut (5o) are cutoff values chosen according to the contours of AT pk (5o, x; Of), 
such as to include only the N peaks that produce the largest temperature shift. The corre- 
sponding calculation for the expected signal for voids simply reverses the sign of the result. 
The variance of the distribution of temperatures about this expectation value can be esti- 
mated by ((AT 2 ) - (AT) 2 ). 

It is worth pointing out that the sample of structures over which the expectation value 
(3.3) is defined is not chosen to match the actual sample of structures in the observed LRG 
density field. Rather it includes only those structures within the specified volume which 
give the largest contribution to the integral, irrespective of whether the particular structure- 
identification algorithm used would find them or not. In other words, (AT(iV, Op)) represents 
the expectation of the maximum possible observable signal from spherical structures, and is 
expected to be achievable only if the structure- finding algorithm is 100% efficient. 

Note also that the procedure for calculation of {AT(N,0p)) outlined above is slightly 
different to that used previously in ref. [28]. In that paper, the radius of a structure of given 
5o was fixed to be at the mean value of the distribution by setting x to its expectation value 
(x\5q). Here we account for the distribution in radii and restrict the sample to include only 
those structures which produce the largest temperature shift. Thus the present analysis is 
even more conservative than that in ref. [28]. 

For N = 50 and Of = 4° and the survey volume V chosen to match that of the 
SDSS DR6 catalogue, the spherical model gives a maximum expected average ISW signal 
of superstructures to be (AT) = 2.27 ± 0.14 fiK, slightly larger than in [28]. Nevertheless, 
comparing this value with the observed value of AT = 9.6 ± 2.2 /uK, we see that the original 
estimate of a discrepancy between theory and observation at a level higher than 3a remains. 
Indeed, the spherical model prediction for the maximum ISW signal remains at the level of 
the observational noise. 

3.2 Limitations of the model 

Although the spherical model described above aims to provide the most optimistic estimate 
for the size of the stacked ISW signal from superstructures, it necessarily rests on certain 
simplifying assumptions. If one or more of these assumptions were to have a large effect on 
the result, it may be that the calculation underestimates the maximum expected ISW signal. 
In such a case the discrepancy between the ACDM model and the observation of ref. [24] 
might be mitigated. 

The first assumption of the model is to place the centres of all structures at the same 
comoving distance from the observer without accounting for any effects of the redshift distri- 
bution. In particular this does not allow for overlap of structures at different redshifts along 
the line of sight. It was originally argued [28] that since overlapping structures could both 
increase or decrease the total signal, depending on the sign of the density fluctuation, the 
net effect of this simplification would be small. 
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Another simplification is that the model uses the linear approximation for growth of 
structures. Although on such large scales and low redshifts one expects this to be valid in 
the mean, with non-linear dynamics providing only a small correction [29], for the specific 
population of the most extreme fluctuations it may not hold. Such non-linear effects would 
be expected to break the symmetry of the predicted signal from voids and peaks. 

The entire calculation is also based on the assumption of Gaussianity of the density 
perturbations. However, non-linear evolution under gravity itself induces a small skewness 
and kurtosis to the density field. Although negligible in the mean, this may be important for 
the extreme tails of the distribution investigated here. Any effect of skewness on the actual 
temperature distribution would also affect voids and peaks differently. 

The assumption of sphericity in particular may lead to an underestimate of the maxi- 
mum possible signal. Although the mean profiles over all structures are spherical, the specific 
population of structures that produce the maximum ISW signal is significantly biased towards 
structures elongated along the line-of-sight direction. Photometric redshift errors in the LRG 
sample cause a significant redshift smearing along the line-of-sight in the reconstructed 3D 
LRG density distribution. As a result of this smearing, the population of structures found by 
structure-finding algorithms will not have isotropic deviations from sphericity. Instead it will 
be somewhat biased towards structures elongated along the line-of-sight direction, which are 
also those that produce the largest ISW temperature shift. This is despite the fact that the 
structure- finding algorithms themselves do not have any preferred directionality. As a result, 
despite the conservative analysis, it is in principle possible that the sample of superstructures 
actually chosen for analysis gives a larger average effect than the model predicts. 

Finally, the spherical model considers all N structures to be statistically independent. 
This is well justified for small values of N; however as N increases this assumption becomes 
less valid. For large values N, the fact that N — 1 extreme structures with specific values of 
5 and x already exist within a given volume will necessarily alter the expected 5 and x for 
the N th structure. 

The net effect of all these assumptions on the spherical model prediction is not clear. 
Indeed the relative importance of each is expected to depend on the redshift interval under 
consideration. What we wish to determine is whether the calculated value (AT(N,0f)) for 
a population of spherical structures can underestimate the true maximum ISW signal that 
might be observed on the sky with a perfect (and lucky) observation. To answer this in the 
next section we compare the model predictions with the maximal stacked signal obtained 
from simulated sky maps of the ISW effect. 

4 Simulated sky maps 

The ISW effect in a ACDM cosmology, including non-linear Rees-Sciama (RS) contributions 
[39], has been studied using large A-body simulations in [29]. This study employed an N- 
body simulation of 2200 3 particles in a box of side 1 /i _1 Gpc, for a ACDM model with 
parameters fi A = 0.74, Sl m = 0.26, fi b = 0.044, a s = 0.8 and H = 71.5 km s" 1 Mpc" 1 . 

For our purposes here a brief summary of the method used to create the simulations will 
suffice. The initial density conditions were set at redshift z = 49 and then evolved through 
50 snapshots to redshift z = 0, with the intervals between neighbouring simulation outputs 
corresponding to a radial comoving distance of ~ 100 /i _1 Mpc. At each of these redshifts, 
the simulation output was used to construct <j> according to eq. (2.3). By interpolating 
between these $ values, photon paths were then traced back through the simulation in order 
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Figure 1. Sky maps of the ISW effect due to the matter density distribution in iV-body simulations 
[29]. The method of generation of the maps is summarised in the text. Map A (left panel) is for the 
rcdshift range 0.17 < z < 0.57, and Map B (right panel) for 0.57 < z < 1.08. The corresponding 
comoving radial distance ranges are indicated in the figure. Note the asymmetry in the definition of 
the colour scale. 



500-1500 Mpc/h, filtered 1500-2500 Mpc/h, filtered 
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Figure 2. The same maps as in figure 1, but with the compensated top-hat filter W(6; Bp) applied. 
The filter radius is Op — 4° in this example. 



to calculate the cumulative AT effect via eq. (2.1) over depths up to the simulation depth of 
1 /i -1 Gpe. When the photon path exits the simulation box, it is mapped back to a location 
within the box using periodic boundary conditions. Finally the AT map thus created is 
visualised in spherical coordinates using HEALPix [40]. 

In this way, full-sky maps of the ISW temperature shift including RS contributions are 
produced. Two of these maps, for the redshift intervals 0.17 < z < 0.57 (corresponding 
to radial comoving distances of 500 — 1500/i _1 Mpc) and 0.57 < z < 1.08 (corresponding 
to 1500 — 2500 Mpc), are shown in figure l. 3 We refer to these as Map A and Map B 
respectively. Although some overlap is present, neither redshift range corresponds exactly 
to that of the SDSS LRG sample described in section 2.2, 0.4 < z < 0.75. Therefore results 
obtained from these maps cannot be directly compared to observed values. However, they 
can be compared to the theoretical predictions for the same redshift intervals made using the 
spherical model of section 3. 

3 We are grateful to Yan-Chuan Cai for providing us with these maps. 
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Figure 3. Magnified section of Map A, illustrating the procedure for selecting the 50 most extreme 
hot and cold pixels within a given window. The large circle defines the window area of ~ 8000 
square degrees. The smaller circular patches within are 4° circles centred on the selected pixels. Red 
circles indicate hot spots and blue circles cold spots. In the left-hand panel, the patches are shown 
superimposed on the unfiltered map, in the right-hand panel on the same area of the filtered map. 
Note that the actual selection is always performed on the filtered map. 



As can be seen in figure 1, the magnitude of the temperature fluctuations are larger 
for Map A, which covers smaller redshifts. This is both because the effect of A on the time 
evolution of <1? is larger at later times, and also because structures are more pronounced. It 
can also be seen that Map A has more power on larger scales. This is expected in any case 
due to the hierarchical growth of density perturbations; however, due to the finite box size 
and the periodic boundary conditions both maps also lack power and show non-Gaussianities 
on angular scales larger than the angle subtended by the 1 /i _1 Gpc box size [29]. This scale 
is naturally smaller at larger comoving distances. 

4.1 Method of analysis 

We wish to use these two sky maps of the ISW effect to test the spherical model prediction 
for the maximum possible stacked signal of superstructures. To do this, we simply work 
backwards, by applying the methodology of ref. [24] to the hottest and coldest spots that can 
be found in the given map, without any reference to actual structures observed in any galaxy 
survey. Clearly the value of AT thus obtained is not a realistic estimate of the expectation 
from an actual observation, but it does provide a robust upper bound to the observable 
stacked signal in ACDM. 

First, we apply the compensated top-hat filter of eq. (2.7) to each pixel of the original 
maps. For a choice of filter radius Op = 4° the resultant filtered maps can be seen in figure 2 
and compared with the equivalent unfiltered maps in figure 1. 

We then select a circular window on the filtered map of angular size ~ 8000 square 
degrees. The size of this window is chosen roughly to match the size of the SDSS DR9 
footprint on the sky, and is somewhat larger than the footprint of the DR6 survey used in 
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Figure 4. An example histogram of the average AT values obtained for 10000 random selections 
of 50 4- 50 pixels on the filtered WMAP ILC map within a circular window of area ~ 8000 square 
degrees. The filter radius is dp — 4° and the average calculated as described in the text. The best-fit 
Gaussian distribution to the histogram is also plotted. 



[24]. 4 Within this window, we select the N hottest and N coldest pixels. 

In order to ensure that not all pixels are chosen from the same structure on the map, we 
apply the following algorithm to limit overlap. For hot spots, we start by selecting the hottest 
pixel within the window after filtering on scale Of- We then remove the circular region of 
radius 9f about this direction from consideration and choose the next hottest pixel, before 
repeating the procedure to obtain N hot spots. The same procedure applied to the coldest 
pixels gives N cold spots. This ensures that the selected pixels are always separated by at 
least Of, but some overlap of the circular regions over which the top-hat filter contributes 
to the pixel temperature is still allowed. Figure 3 illustrates the result of one such selection 
process. Averaging the temperatures of these pixels gives the quantities ATh t, AT co id and 
AT = (AThot — AT co id)/2 as functions of N and 6p for the given location of the circular 
window on the sky. 

To account for the variation over the whole map, we repeat the measurement for 10 
different locations of the window and calculate the average AT(N,0f) and the variance 
<j\ t (N,0f) over these locations for hot spots, cold spots and the total separately. The 
window locations are chosen to ensure coverage of the entire map. As a result there is some 
overlap between different windows and the AT values obtained from each are not completely 
independent. This overlap is increased if more windows are used, leading to an artificial 
underestimate of the variance. However, we checked that the mean AT(N,0f) remains 
stable. It should also be noted that each map is constructed from only one 1 /i _3 Gpc 3 
realisation of the density field, which needs to be used more than once to generate a map of 
the entire sky. Nevertheless, the mean AT(N, Of) obtained is an unbiased estimator of the 
true ensemble average value. 

We emphasize again that AT(N, Of) obtained from the maps in this fashion is a measure 

4 Of course, the circular shape of the window does not correspond to the actual SDSS footprint, but the 
difference due to the geometry is small. We repeated the analysis with other shapes of the window and found 
no significant change to the result. 
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Figure 5. The solid curves show the variation of the maximum AThot (red) and AT co id (blue) as 
a function of filter radius Op as determined from the simulated maps with N = 50. The approximate 
error bars are obtained from the variation due to the choice of location of the window area on the 
maps. The dashed curves show the corresponding theoretical maximum values (AT C ) and (AT V ) 
obtained from the spherical model. The left panel is for Map A and the right panel for Map B. 



of the maximum stacked temperature signal that could be in principle be observed in a 
ACDM cosmology by an ideal experiment in which the hot and cold spots were chosen based 
on an a posteriori knowledge of the actual ISW temperature distribution. It should not be 
regarded as representative of the expected signal from any realistic observation. 

In order to assess the significance of the average maximum signal thus found, we need to 
estimate the noise introduced in the measurement by primary CMB fluctuations. To do this 
we use the WMAP 7-year Internal Linear Combination map 5 and apply the compensated 
top-hat filter of radius Of- We then choose 2N random pixels from within the same circular 
window as before, placed at a location chosen to maximize the overlap with the actual SDSS 
DR6 footprint, using the procedure described above to limit the amount of overlap. We add 
the temperatures of the first N pixels, subtract the temperatures of the next N pixels, and 
divide the result by 2N to obtain one random realization of AT. We repeat this procedure 
10, 000 times and fit a Gaussian distribution function to the resulting histogram of AT. The 
normalized histogram and fitted distribution function for N = 50 and Op = 4° are shown in 
figure 4. The standard deviation of the fitted distribution provides the uncorrelated noise 
on the stacked ISW temperature measurement. The noise values we obtain agree well with 
those reported in the original measurement [24]. 

5 Results and model comparison 

Figure 5 shows the variation of AT(N,Qp) as a function of filter radius Of for N = 50 
obtained from the simulated maps for hot and cold spots separately. This represents the 
maximum possible ISW temperature signal in a ACDM universe. Also plotted is the theoret- 
ical prediction for the spherical model described in section 3. The dependence of AT(N, Of) 
on N for Of = 4° is shown in figure 6. 



5 Available to download from http://lambda.gsfc.nasa.gov/product/map/dr4/ilc_map_info.cfm. 
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Figure 6. Same as figure 5 except showing the variation of the maximum AT as a function of N, 
for fixed filter radius 6f — 4°. 



Comparing the spherical model predictions and the upper limits AThot/coid obtained 
from Map A, which covers comoving distances in the range 500 — 1500 fa -1 Mpc, we find 
that the spherical model somewhat underestimates the maximum possible signal that is in 
principle measurable by an ideal experiment. The values of AThot an d AT co id obtained from 
the maps are almost symmetrical, as would be expected for structures in the linear regime. 
This indicates that neither the non-linear evolution of density profiles nor the induced non- 
Gaussianity of the density distribution have any significant effect on the stacked ISW signal 
on the relevant scales. 

Instead the reason for the difference in Map A is probably the assumption of sphericity. 
It is likely that the structures in the ./V-body simulation responsible for the very largest 
temperature shifts in this map are better described as elongated filamentary structures that 
happen to be aligned along the line of sight. By approximating them as spheres, the spherical 
model underestimates their ISW temperature shift. Since such long filaments that happen 
to be aligned along the line of sight are rare, their influence on the average AT is relatively 
smaller when larger N populations are considered, hence the amount by which the model 
prediction underestimates the actual maximal signal in Map A decreases at large N. The 
effect of neglecting correlations between extreme structures would also result in an artificial 
enhancement in the spherical model's prediction for large N. 

It should be emphasized again that the spherical model prediction is for the expectation 
of the ISW temperature signal associated with peaks in the density fluctuation field, whereas 
the AThot/coid values are associated with peaks in the potential field and may not be directly 
observable. Therefore AT hot / cold represent an upper bound on the maximum possible stacked 
ISW signal. Although for Map A the spherical model underestimates this maximum possible 
signal for the reasons stated above, the numerical difference is still small for all ./V and 6p. 
In particular, at N = 50 and Op = 4° the difference is < 0.5 /xK. 

On the other hand for Map B, which covers comoving distances in the range 1500 — 
2500 /i _1 Mpc, the spherical model prediction for the maximum possible signal is slightly 
larger than the map values. At these distances, relatively narrow filamentary structures 
subtend a much smaller angle on the sky and so produce a smaller effective AT signal after 
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application of the filter of width Op. The structures that give the largest effective AT will be 
much wider, and therefore the spherical approximation is less important. Instead, effects such 
as the correlation between extreme structures mentioned above appear to be more relevant. 
As the spherical model overestimates the maximum possible signal at these high redshifts it 
can be used to derive a conservative upper bound on the observable signal. 

Another quantity of interest is the ratio of the maximum possible stacked signal obtained 
from the map analysis to the intrinsic observational noise in the measurement due to primary 
CMB anisotropies. In figures 7 and 8 we plot the behaviour of the signal-to-noise ratio as 
a function of N and Op. In this case while varying 9p the pixel selection is always based 
on the map filtered with the 4° filter in order to keep the same sample of structures. The 
signal-to- noise ratio does not show a peak at N = 50 or at Op = 4° for either of the maps 
analyzed, in contrast to the observation in ref. [24]. It is also clear that even though by 
construction the map analysis selects precisely the regions contributing the largest filtered 
signal, the maximum signal-to-noise ratio in a ACDM cosmology with Gaussian perturbations 
is always less than 1.5. A realistic observation would not even be able to achieve this value. 
We conclude that the observed 4.4<r significance of the signal is totally at odds with the 
concordance cosmology. 

6 Conclusions 

The results of the previous section show that at high redshifts, the spherical model calculation 
for the maximum possible stacked ISW signal that could be observed is an overestimate of 
the actual value obtained from the simulated map as described. At lower redshifts the a 
posteriori map analysis gives a larger maximum signal than the spherical model calculation, 
though the quantitative differences are small. 

The symmetry between hot and cold spots in the simulated maps shows that the linear 
treatment in the spherical model is valid and that the non-Gaussianity induced by gravita- 
tional evolution does not produce a significant effect. 

Instead at low redshifts the spherical model underestimates the AT from the simulated 
map primarily due to the assumption of sphericity. Whether the actual maximum possible 
stacked ISW signal in a realistic observation is closer to the spherical model prediction or 
the value from the map analysis depends on the properties of the galaxy catalogue and the 
strategy that is used to identify the superstructures. For a catalogue with spectroscopically 
determined redshifts, the redshift smearing problem described in section 3.2 is less severe. 
Therefore a structure-finding algorithm that does not have a directional preference will se- 
lect a population of structures that, although not spherical, will have uniformly distributed 
orientations and whose mean effect can therefore be well described by the spherical model. 
On the other hand when only photometric redshifts are available, the redshift smearing in- 
duced will naturally bias the population of structures selected towards including more long 
structures aligned along the observer's line of sight. In such a case despite its conservative 
assumptions, the spherical model underestimates the maximum stacked signal that can in 
principle be observed. 6 

Of course the closeness of the two predictions at low redshifts mean that this is a 
moot point. The values of (AT) and AT obtained from Map A for the redshift range 
0.17 < z < 0.57 differ by less than 0.5 (jlK at N = 50 and 9 F = 4°. For Map B, over 

6 We note in passing that in the presence of such redshift errors, searching for structures in the reconstructed 
3D density field is likely not the optimal strategy in any case. 
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Figure 7. The variation of the maximum achievable signal-to-noise ratio as a function of filter 
radius Op, for fixed N = 50. The maximum signal is calculated by subtracting the average of the 
cold spots from the average of hot spots and dividing by 2N, and the errors are added in quadrature. 
The noise is calculated from the width of the Gaussian fit to the distribution of AT values found for 
a random location of patches on the WMAP ILC map, as described in detail in the text. 
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Figure 8. Same as figure 7, except showing the maximum achievable signal-to-noise ratio as a 
function of N, for fixed 9 F =4°. 



the redshift range 0.57 < z < 1.08, the difference is even smaller and the spherical model 
overestimates the obtainable signal. 

We therefore conclude that the spherical model prediction (AT) = 2.27±0.14 fiK using 
the actual SDSS DR6 survey volume and redshift distribution can meaningfully be compared 
with the observed value of AT = 9.6±2.2 fiK for the combined signal of voids and clusters [24] . 
The observed signal is easily > 3a larger than even the most optimistic expectation in the 
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standard ACDM cosmological model with Gaussian primordial perturbations. This confirms 
the conclusion reached in ref. [28] that the observed correlation cannot be due to the linear 
ISW effect in such a cosmology. 

The question that then arises is, what does account for this large correlation between 
structures in the LRG sample and CMB temperature fluctuations? The immediate assump- 
tion is that it is a spurious signal caused by some as yet undetected systematic or foreground 
contamination in the observation. It is however difficult to explain how this might be the 
case. Contamination in the CMB maps from extragalactic radio sources or due to Sunyaev- 
Zel'dovich (SZ) effect is potentially possible. However, the original observation method used 
CMB maps with the conservative KQ75 mask applied and foreground templates subtracted 
to avoid contamination from the Galaxy and point sources. The authors also checked for, but 
did not find, any frequency dependence of the signal over the Q, V and W band CMB maps 
in the absence of foreground subtraction [24]. This suggests that foreground contamination 
is not an issue. 

The SZ effect would be expected to reduce the signal strength for the hot spots associated 
with overdensities, though the few clusters present in voids might enhance the signal seen 
for cold spots. The characteristic scale of this effect is however small compared to the filter 
radius used here, and so the net effect would be small. Nevertheless, it would be a worthwhile 
exercise to repeat the same observational method using a mask from the improved SZ point 
source catalogue available from Planck. 

If foreground contamination is not important and no other systematic effect can be 
found, it is still possible that the observation is simply due to an unfortunate correlation 
between primary CMB fluctuations and the directions along which superstructures were 
identified in the LRG catalogue. Such a correlation is unlikely at a > 3a level given the 
size of the ISW signal expected in ACDM (if selecting CMB patches at random without 
any reference to large-scale structure it is > 4.4a unlikely), but even such rare fluctuations 
do of course sometimes occur. However, if the signal is caused simply by such an unlucky 
fluctuation, it should disappear if the same observation is repeated with another catalogue 
of galaxies or other tracers of the matter density field that covers a different region of the 
sky. Efforts to test for such a correlation are somewhat complicated by the fact that in the 
absence of accurate redshift information for the given catalogue, the method used here to 
locate structures in the 3D distribution is not optimal and other techniques must be devised. 
This is work in progress. 

If such further analyses should confirm the existence of the stacked ISW signal then 
some cosmological explanation is required. This would be a strong indication of the need 
for some new physics beyond the standard ACDM model. We have previously speculated 
as to what such new physics might be [28]. Taken at face value, the observation indicates 
that the rarest and most extreme fluctuations in the matter density contrast field 5 are in 
fact more numerous than expected. This means that the extreme tails of the probability 
distribution function for 5 are modified from the Gaussian expectation on linear scales. Such 
a modification could arise from a primordial non-Gaussian kurtosis in the matter density field, 
though the magnitude of the kurtosis required to explain the results is still to be calculated. 
Alternatively, the explanation could lie in the late-time growth of structures in an alternative 
to ACDM, such as modified gravity or f{R) models. 

Such a modification of the Gaussian distribution function of 5 would also affect the 
full cross-correlation amplitude for the ISW effect measured in other studies, although such 
observations are not specifically designed to test the tail of the distribution so might be 
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less sensitive to deviations. We note again that mild excesses of the full cross-correlation 
amplitude have sometimes been found when using the same LRG catalogues as tracers [19- 



We conclude by noting that the stacked ISW signal of rare superstructures is one of the 
most statistically significant anomalies that has yet been found in the standard cosmological 
model. The ISW effect is a linear order effect and in this work we have demonstrated 
that the theoretical calculation of the expectation value on relevant scales is not affected by 
complexities of non- linear dynamics. It is also independent of any assumptions about the 
precise nature of the bias relationship between the galaxy distribution and the matter density 
field. This may be taken as further reason to believe in the robustness of this anomaly. 
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Note added: Shortly after submission of this paper, an independent study of the same issue 
was uploaded to the arXiv [41]. This study uses direct A-body simulations of the relevant 
redshift range and finds the same signal size for ACDM as presented here. The authors also 
re-perform the stacking analysis and confirm that the signal seen in [24] is robust to rotation 
tests of the CMB maps, while finding no evidence for foreground contamination. 

A Statistics of Gaussian random fields 

In this appendix, we provide definitions for some of the quantities and functions used in 
section 3. A full derivation of these expressions from a treatment of the statistics of Gaussian 
random fields is provided in ref. [34]; here we limit ourselves to the details required for 
completeness of this paper. 

Consider a homogeneous and isotropic Gaussian random field S(r), with two-point cor- 
relation function £(|ri — 1"2 1 ) = (<5(ri)<5(r2)) and power spectrum P(k) defined as the Fourier 
transform of £(r). The field is smoothed on length scale R{ using a Gaussian filter, such that 
in Fourier space the smoothed field and its power spectrum are 



22]. 



5f(k) = exp (-R 2 k 2 /2) 6(k) , P t {k) = exp(-R 2 k 2 )P{k) . 



(A.l) 



The spectral moments of this field, u,-, are defined as 




(A.2) 



and the corresponding spectral parameters are 



7 = 




(A.3) 
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Throughout the paper for the purposes of numerical evaluation we take a smoothing scale of 
Rf = 20 h^Mpc. 

The expression for the differential comoving number density of peaks of the smoothed 
field as a function of v and x can then be expressed in terms of these spectral parameters as 
in eq. (3.1). The function f(x) appearing in this expression can be written in closed form as 
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The mean spherical profile about a peak with given v and x can be written as 



5 pk (r;z/,x) 
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7 2 ^ + ^VV 



(A.4) 



(A.5) 



where ip(r) = C( r )/?(0)- Writing £(r) in terms of its Fourier transform, using £(0) = <rg and 
rearranging terms leads to eq. (3.2). 
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